09 Resumen regresión dinámica

Author

Eddie Aguilar

library("easypackages")
packages("tidyverse","fpp3", "tsibble", "feasts","fable", "patchwork")
library(plotly)

Modelos de regresión dinámica

1. Introducción

Los modelos de suavización exponencial y ARIMA no pueden incluir información externa a la serie. Por el otro lado, regresión lineal sí toma en cuenta otras variables, pero siempre asume que el error es ruido blanco. Por lo tanto, podemosunir estos dos modelos y resultará en un modelo de regresión dinámica, igual que al lineal pero el error es modelado con ARIMA.

Hay que tomar en cuenta que en regresión dinámica, es necesario que todas las variables sean estacionarias.

2. Regresión con errores ARIMA en R con fable

Al agregar pdq(d=1), R aplicará las primeras diferencias a todas las variables, si no se agrega los valores de pdq, R los buscará automáticamente.

ARIMA(y ~ x + pdq(1,1,0))
<ARIMA model definition>

Ejemplo

Analizar y pronosticar los cambios en el consumo personal a través de el ingreso disponible.

us_change <- read_csv("https://otexts.com/fpp3/extrafiles/us_change.csv") %>%
  mutate(Time = yearquarter(Time)) %>%
  as_tsibble(index = Time)
Rows: 187 Columns: 6
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl  (5): Consumption, Income, Production, Savings, Unemployment
date (1): Time

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
us_change %>%
  gather("var", "value", Consumption, Income) %>%
  ggplot(aes(x = Time, y = value)) +
  geom_line() +
  facet_grid(vars(var), scales = "free_y") +
  xlab("Year") + ylab(NULL) +
  ggtitle("Quarterly changes in US consumption and personal income")

us_change %>% 
  ggplot(aes(x = Income, y = Consumption))+
  geom_point()

pdq automáticos:

fit <- us_change %>%
  model(ARIMA(Consumption ~ Income))
report(fit)
Series: Consumption 
Model: LM w/ ARIMA(1,0,2) errors 

Coefficients:
         ar1      ma1     ma2  Income  intercept
      0.6922  -0.5758  0.1984  0.2028     0.5990
s.e.  0.1159   0.1301  0.0756  0.0461     0.0884

sigma^2 estimated as 0.3219:  log likelihood=-156.95
AIC=325.91   AICc=326.37   BIC=345.29

Viendo los errores de la regresión y de ARIMA con residuals(), los podemos seperar con type = “regression” y type = “innovations” para los errores ARIMA.

bind_rows(
  `Regression Errors` = as_tibble(residuals(fit, type="regression")),
  `ARIMA Errors` = as_tibble(residuals(fit, type="innovation")),
  .id = "type"
) %>%
  ggplot(aes(x = Time, y = .resid)) +
  geom_line() +
  facet_grid(vars(type), scales = "free_y") +
  xlab("Year") + ylab(NULL)

Asegurarnos que los errores ARIMA sean ruido blanco:

fit %>% gg_tsresiduals()

Ljung-Box:

augment(fit) %>%
  features(.resid, ljung_box, dof = 5, lag = 8)
# A tibble: 1 × 3
  .model                      lb_stat lb_pvalue
  <chr>                         <dbl>     <dbl>
1 ARIMA(Consumption ~ Income)    5.89     0.117

3. Pronóstico

Para llevar a cabo pronósticos de modelos de regresión con errores ARIMA, se necesita realizar el pronóstico de:

  • La parte de la regresión
  • La parte de los erroes ARIMA

y combinar los resultados.

Necesitamos pronósticos de las variables independientes o predictoras para poder pronosticar nuestra variable de interés. Si las predictoras no son conocidas en el futuro, tenemos que modelarlas por separado, o asumir valores futuros para cada una.

Ejemplo:

Regresando al consumo personal a partir del ingreso, podemos suponer que los cambios porcentuales en el ingreso serán iguales a el cambio promedio porcentual de los últimos 40 años:

us_change_future <- new_data(us_change, 8) %>% mutate(Income = mean(us_change$Income))
forecast(fit, new_data = us_change_future) %>%
  autoplot(slice_tail(us_change, n = 80)) + 
  xlab("Year") +
  ylab("Percentage change") + 
  ggtitle("Pronóstico de regresión con errores ARIMA")

Comparando con ARIMA:

fit_prev <- us_change %>%
  model(ARIMA(Consumption ~ PDQ(0,0,0)))

fit_prev %>% forecast(h=10) %>% 
  autoplot(slice(us_change, (n()-80):n())) + 
  ggtitle("Pronóstico con modelo ARIMA")

En este nuevo modelo, podemos capturar más informción, por lo tanto, los intervalos de predicción se reducen.

Hay que tomar en cuenta que los intervalos de predicción no toman en cuenta la incertidumbre de las predictoras.

Ejemplo:

Demanda de energía a parir de temperatura.

vic_elec_daily <- vic_elec %>%
  filter(year(Time) == 2014) %>%
  index_by(Date = date(Time)) %>%
  summarise(
    Demand = sum(Demand)/1e3,
    Temperature = max(Temperature),
    Holiday = any(Holiday)
  ) %>%
  mutate(Day_Type = case_when(
    Holiday ~ "Holiday",
    wday(Date) %in% 2:6 ~ "Weekday",
    TRUE ~ "Weekend"
  ))

vic_elec_daily %>%
  ggplot(aes(x=Temperature, y=Demand, colour=Day_Type)) +
    geom_point() +
    ylab("Electricity demand (GW)") +
    xlab("Maximum daily temperature")

Como vemos, es una función cuadratica y la energia demandada sí depende del día y la temperatura.

Gráficas de tiempo:

vic_elec_daily %>% 
  pivot_longer(cols = c(Demand, Temperature), names_to = "vars",
               values_to = "value") %>% 
  ggplot(aes(x = Date, y = value)) + 
  geom_line() + 
  facet_wrap(~ vars, ncol = 1, strip.position = "right", scales = "free")

Al ser una relación cuadrática, ajustaremos el modelo cuadrático con errores ARIMA, y agregamos una variable para indicar si el día fue hábil o no.

fit <- vic_elec_daily %>%
  model(ARIMA(Demand ~ Temperature + I(Temperature^2) + (Day_Type=="Weekday")))

report(fit)
Series: Demand 
Model: LM w/ ARIMA(2,1,2)(2,0,0)[7] errors 

Coefficients:
          ar1     ar2      ma1      ma2    sar1    sar2  Temperature
      -0.1093  0.7226  -0.0182  -0.9381  0.1958  0.4175      -7.6135
s.e.   0.0779  0.0739   0.0494   0.0493  0.0525  0.0570       0.4482
      I(Temperature^2)  Day_Type == "Weekday"TRUE
                0.1810                    30.4040
s.e.            0.0085                     1.3254

sigma^2 estimated as 44.91:  log likelihood=-1206.11
AIC=2432.21   AICc=2432.84   BIC=2471.18
fit %>% gg_tsresiduals()

augment(fit) %>%
  features(.resid, ljung_box, dof = 8, lag = 14)
# A tibble: 1 × 3
  .model                                                         lb_stat lb_pv…¹
  <chr>                                                            <dbl>   <dbl>
1 "ARIMA(Demand ~ Temperature + I(Temperature^2) + (Day_Type ==…    28.4 7.91e-5
# … with abbreviated variable name ¹​lb_pvalue

Al momento de pronosticar, necesitamos datos respecto a la tempreatura futura y podemos hacer dos cosas:

  • Conseguir los datos de algún pronóstico meteorológico y meterlos al modelo.
  • Podemos pronosticar bajo un escenario. Asumiremos que la temperatura máxima se mantendrá constante a 26 grados.
vic_elec_future <- new_data(vic_elec_daily, 14) %>%
  mutate(
    Temperature = 26,
    Holiday = c(TRUE, rep(FALSE, 13)),
    Day_Type = case_when(
      Holiday ~ "Holiday",
      wday(Date) %in% 2:6 ~ "Weekday",
      TRUE ~ "Weekend"
    )
  )
forecast(fit, vic_elec_future) %>%
  autoplot(vic_elec_daily) + ylab("Electricity demand (GW)")

vic_elec_future <- new_data(vic_elec_daily, 14) %>%
  mutate(
    Temperature = c(seq(31,34),33,rep(34,3),rep(33,6)),
    Holiday = c(TRUE, rep(FALSE, 13)),
    Day_Type = case_when(
      Holiday ~ "Holiday",
      wday(Date) %in% 2:6 ~ "Weekday",
      TRUE ~ "Weekend"
    )
  )
vic_elec_future
# A tsibble: 14 x 4 [1D]
   Date       Temperature Holiday Day_Type
   <date>           <dbl> <lgl>   <chr>   
 1 2015-01-01          31 TRUE    Holiday 
 2 2015-01-02          32 FALSE   Weekday 
 3 2015-01-03          33 FALSE   Weekend 
 4 2015-01-04          34 FALSE   Weekend 
 5 2015-01-05          33 FALSE   Weekday 
 6 2015-01-06          34 FALSE   Weekday 
 7 2015-01-07          34 FALSE   Weekday 
 8 2015-01-08          34 FALSE   Weekday 
 9 2015-01-09          33 FALSE   Weekday 
10 2015-01-10          33 FALSE   Weekend 
11 2015-01-11          33 FALSE   Weekend 
12 2015-01-12          33 FALSE   Weekday 
13 2015-01-13          33 FALSE   Weekday 
14 2015-01-14          33 FALSE   Weekday 
forecast(fit, vic_elec_future) %>%
  autoplot(vic_elec_daily) + ylab("Electricity demand (GW)")

4. Regresión armónica dinámica

Las series de Fourier es un método alternativo de los dummies muy bueno para periodos estacionales muy largos. Ya que ARIMA y ETS son métodos que estiman estacionalidad, pero pierde efectividad entre mayor sea el periodo estacional, y que, las dummies estiman estacionalidad de cualquier tamaño, pero si es un periodo muy grande, se necesitarán muchas variables, el uso de series de Fourier es buena opción.

Una serie de Fourier es una función básica trigonométrica, tiene un seno y un coseno. Al ser una función trigonométrica periódica, simula una estacionalidad y el uso de esta nos puede ayudar a pronosticar una serie estacional.

Como se comentó, una serie de Fourier se consituye de un seno y coseno, y se pueden utilizar tantos pares de senos y cosenos como se desea.

El número de pares es lo que cambiará la suavización o el ajuste del método a los datos de entrenamiento. Entre más pares, mejor ajuste (hay que tomar en cuenta el over-fitting).

El máximo absoluto de pares usados para un buen pronóstico es el periodo estacional \(m / 2\).

Ejemplos de estacionalidades largas:

  • Datos diarios, donde pueden tener estacionalidad anual. (m = 365)
  • Datos semanales (m = 52)
  • Datos por cada media hora con estacionalidad diaria. (m = 48)

Ventajas:

  • Permite cualquier tamaño de estacionalidad.
  • Se pueden simular varios tipos de estacionalidad al mismo tiempo (con distintas frecuencias).
  • La suavización del patrón estacional es controlado.
  • La dinámica de corto plazo puede ser controlada a través de errores ARIMA.

Ejemplo

Gasto en comidas fuera de casa.

aus_cafe <- aus_retail %>%
  filter(
    Industry == "Cafes, restaurants and takeaway food services",
    year(Month) %in% 2004:2018
  ) |> summarise(Turnover = sum(Turnover))
aus_cafe
# A tsibble: 180 x 2 [1M]
      Month Turnover
      <mth>    <dbl>
 1 2004 Jan    1895.
 2 2004 Feb    1766.
 3 2004 Mar    1873.
 4 2004 Apr    1874.
 5 2004 May    1846 
 6 2004 Jun    1758.
 7 2004 Jul    1878.
 8 2004 Aug    1850.
 9 2004 Sep    1922.
10 2004 Oct    1915.
# … with 170 more rows
autoplot(aus_cafe)
Plot variable not specified, automatically selected `.vars = Turnover`

Probamos con varios niveles \(K\) (número de pares):

fit <- aus_cafe %>%
  model(
    `K = 1` = ARIMA(log(Turnover) ~ fourier(K = 1) + PDQ(0,0,0)),
    `K = 2` = ARIMA(log(Turnover) ~ fourier(K = 2) + PDQ(0,0,0)),
    `K = 3` = ARIMA(log(Turnover) ~ fourier(K = 3) + PDQ(0,0,0)),
    `K = 4` = ARIMA(log(Turnover) ~ fourier(K = 4) + PDQ(0,0,0)),
    `K = 5` = ARIMA(log(Turnover) ~ fourier(K = 5) + PDQ(0,0,0)),
    `K = 6` = ARIMA(log(Turnover) ~ fourier(K = 6) + PDQ(0,0,0))
  )
glance(fit)
# A tibble: 6 × 8
  .model   sigma2 log_lik   AIC  AICc   BIC ar_roots  ma_roots 
  <chr>     <dbl>   <dbl> <dbl> <dbl> <dbl> <list>    <list>   
1 K = 1  0.00175     317. -616. -615. -588. <cpl [2]> <cpl [3]>
2 K = 2  0.00107     362. -700. -698. -661. <cpl [5]> <cpl [1]>
3 K = 3  0.000761    394. -763. -761. -725. <cpl [3]> <cpl [1]>
4 K = 4  0.000539    427. -822. -818. -771. <cpl [1]> <cpl [5]>
5 K = 5  0.000317    474. -919. -917. -875. <cpl [2]> <cpl [0]>
6 K = 6  0.000316    474. -920. -918. -875. <cpl [0]> <cpl [1]>
rep_model <- function(modelo) {
  print(strrep("- - ",18))
  print(paste("Modelo", modelo, sep = " "))
  fit %>% 
    select(all_of(modelo)) %>% 
    report()
}

names(fit) %>% 
  map(rep_model)
[1] "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - "
[1] "Modelo K = 1"
Series: Turnover 
Model: LM w/ ARIMA(2,1,3) errors 
Transformation: log(Turnover) 

Coefficients:
          ar1     ar2      ma1      ma2     ma3  fourier(K = 1)C1_12
      -0.2266  0.0471  -0.6371  -0.5810  0.6407               0.0067
s.e.   0.1119  0.1103   0.0817   0.0682  0.0651               0.0018
      fourier(K = 1)S1_12  intercept
                  -0.0358     0.0041
s.e.               0.0018     0.0011

sigma^2 estimated as 0.001747:  log likelihood=317.24
AIC=-616.47   AICc=-615.41   BIC=-587.78
[1] "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - "
[1] "Modelo K = 2"
Series: Turnover 
Model: LM w/ ARIMA(5,1,1) errors 
Transformation: log(Turnover) 

Coefficients:
          ar1      ar2      ar3     ar4     ar5      ma1  fourier(K = 2)C1_12
      -0.7492  -0.6512  -0.0657  0.2450  0.5088  -0.2694               0.0067
s.e.   0.1299   0.1497   0.1650  0.1223  0.0750   0.1442               0.0020
      fourier(K = 2)S1_12  fourier(K = 2)C2_12  fourier(K = 2)S2_12  intercept
                  -0.0363              -0.0044              -0.0217     0.0039
s.e.               0.0020               0.0015               0.0015     0.0010

sigma^2 estimated as 0.001073:  log likelihood=361.85
AIC=-699.71   AICc=-697.83   BIC=-661.46
[1] "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - "
[1] "Modelo K = 3"
Series: Turnover 
Model: LM w/ ARIMA(3,1,1) errors 
Transformation: log(Turnover) 

Coefficients:
          ar1     ar2     ar3      ma1  fourier(K = 3)C1_12
      -0.4021  0.2096  0.4883  -0.6694               0.0067
s.e.   0.2453  0.2633  0.1475   0.2531               0.0022
      fourier(K = 3)S1_12  fourier(K = 3)C2_12  fourier(K = 3)S2_12
                  -0.0363              -0.0045              -0.0215
s.e.               0.0022               0.0014               0.0014
      fourier(K = 3)C3_12  fourier(K = 3)S3_12  intercept
                  -0.0071              -0.0355     0.0039
s.e.               0.0016               0.0016     0.0009

sigma^2 estimated as 0.0007609:  log likelihood=393.61
AIC=-763.21   AICc=-761.33   BIC=-724.96
[1] "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - "
[1] "Modelo K = 4"
Series: Turnover 
Model: LM w/ ARIMA(1,1,5) errors 
Transformation: log(Turnover) 

Coefficients:
          ar1      ma1      ma2     ma3      ma4     ma5  fourier(K = 4)C1_12
      -0.7133  -0.0944  -0.0834  0.1182  -0.2059  0.3678               0.0067
s.e.   0.0806   0.1011   0.0720  0.0927   0.0642  0.0671               0.0018
      fourier(K = 4)S1_12  fourier(K = 4)C2_12  fourier(K = 4)S2_12
                  -0.0364              -0.0043              -0.0215
s.e.               0.0018               0.0019               0.0019
      fourier(K = 4)C3_12  fourier(K = 4)S3_12  fourier(K = 4)C4_12
                  -0.0069              -0.0355               0.0034
s.e.               0.0012               0.0012               0.0019
      fourier(K = 4)S4_12  intercept
                  -0.0203     0.0039
s.e.               0.0019     0.0011

sigma^2 estimated as 0.0005386:  log likelihood=426.78
AIC=-821.57   AICc=-818.21   BIC=-770.57
[1] "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - "
[1] "Modelo K = 5"
Series: Turnover 
Model: LM w/ ARIMA(2,1,0) errors 
Transformation: log(Turnover) 

Coefficients:
          ar1      ar2  fourier(K = 5)C1_12  fourier(K = 5)S1_12
      -0.4448  -0.1193               0.0068              -0.0364
s.e.   0.0741   0.0747               0.0024               0.0024
      fourier(K = 5)C2_12  fourier(K = 5)S2_12  fourier(K = 5)C3_12
                  -0.0044              -0.0216              -0.0070
s.e.               0.0014               0.0014               0.0013
      fourier(K = 5)S3_12  fourier(K = 5)C4_12  fourier(K = 5)S4_12
                  -0.0355               0.0035              -0.0202
s.e.               0.0013               0.0014               0.0014
      fourier(K = 5)C5_12  fourier(K = 5)S5_12  intercept
                  -0.0069              -0.0249     0.0039
s.e.               0.0014               0.0014     0.0008

sigma^2 estimated as 0.0003173:  log likelihood=473.73
AIC=-919.47   AICc=-916.91   BIC=-874.85
[1] "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - "
[1] "Modelo K = 6"
Series: Turnover 
Model: LM w/ ARIMA(0,1,1) errors 
Transformation: log(Turnover) 

Coefficients:
          ma1  fourier(K = 6)C1_12  fourier(K = 6)S1_12  fourier(K = 6)C2_12
      -0.3954               0.0068              -0.0363              -0.0044
s.e.   0.0619               0.0024               0.0024               0.0016
      fourier(K = 6)S2_12  fourier(K = 6)C3_12  fourier(K = 6)S3_12
                  -0.0215              -0.0070              -0.0355
s.e.               0.0016               0.0014               0.0014
      fourier(K = 6)C4_12  fourier(K = 6)S4_12  fourier(K = 6)C5_12
                   0.0035              -0.0202              -0.0069
s.e.               0.0013               0.0013               0.0013
      fourier(K = 6)S5_12  fourier(K = 6)C6_12  intercept
                  -0.0249               0.0014     0.0039
s.e.               0.0013               0.0009     0.0008

sigma^2 estimated as 0.0003163:  log likelihood=474.03
AIC=-920.06   AICc=-917.5   BIC=-875.44
[[1]]
# A mable: 1 x 1
                      `K = 1`
                      <model>
1 <LM w/ ARIMA(2,1,3) errors>

[[2]]
# A mable: 1 x 1
                      `K = 2`
                      <model>
1 <LM w/ ARIMA(5,1,1) errors>

[[3]]
# A mable: 1 x 1
                      `K = 3`
                      <model>
1 <LM w/ ARIMA(3,1,1) errors>

[[4]]
# A mable: 1 x 1
                      `K = 4`
                      <model>
1 <LM w/ ARIMA(1,1,5) errors>

[[5]]
# A mable: 1 x 1
                      `K = 5`
                      <model>
1 <LM w/ ARIMA(2,1,0) errors>

[[6]]
# A mable: 1 x 1
                      `K = 6`
                      <model>
1 <LM w/ ARIMA(0,1,1) errors>
p <- fit %>%
  forecast(h = "2 years") %>%
  autoplot(aus_cafe) +
  facet_wrap(vars(.model), ncol = 2) +
  guides(colour = FALSE) +
  geom_label(
    aes(x = yearmonth("2007 Jan"), y = 4250, label = paste0("AICc = ", format(AICc))),
    data = glance(fit)
  )
Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
of ggplot2 3.3.4.
p

p <- fit %>% 
  select(`K = 3`) %>% 
  forecast(h = "2 years") %>%
  autoplot(aus_cafe) +
  facet_wrap(vars(.model), ncol = 2) +
  guides(colour = FALSE)
p

Ejemplo 2

vic_elec
# A tsibble: 52,608 x 5 [30m] <Australia/Melbourne>
   Time                Demand Temperature Date       Holiday
   <dttm>               <dbl>       <dbl> <date>     <lgl>  
 1 2012-01-01 00:00:00  4383.        21.4 2012-01-01 TRUE   
 2 2012-01-01 00:30:00  4263.        21.0 2012-01-01 TRUE   
 3 2012-01-01 01:00:00  4049.        20.7 2012-01-01 TRUE   
 4 2012-01-01 01:30:00  3878.        20.6 2012-01-01 TRUE   
 5 2012-01-01 02:00:00  4036.        20.4 2012-01-01 TRUE   
 6 2012-01-01 02:30:00  3866.        20.2 2012-01-01 TRUE   
 7 2012-01-01 03:00:00  3694.        20.1 2012-01-01 TRUE   
 8 2012-01-01 03:30:00  3562.        19.6 2012-01-01 TRUE   
 9 2012-01-01 04:00:00  3433.        19.1 2012-01-01 TRUE   
10 2012-01-01 04:30:00  3359.        19.0 2012-01-01 TRUE   
# … with 52,598 more rows

Usando datos cada hora:

(vic_elec_hour <- index_by(vic_elec, time = ~ floor_date(., unit="hour")) |> 
  summarise(demand = sum(Demand),
         mean_temp = mean(Temperature),
         max_temp = max(Temperature),
         holiday = any(Holiday)) |> 
    mutate(
      day_type = case_when(
        holiday ~ "holiday",
        wday(time) %in% 2:6 ~ "weekday",
        TRUE ~ "weekend"
      )))
# A tsibble: 26,304 x 6 [1h] <Australia/Melbourne>
   time                demand mean_temp max_temp holiday day_type
   <dttm>               <dbl>     <dbl>    <dbl> <lgl>   <chr>   
 1 2012-01-01 00:00:00  8646.      21.2     21.4 TRUE    holiday 
 2 2012-01-01 01:00:00  7927.      20.6     20.7 TRUE    holiday 
 3 2012-01-01 02:00:00  7902.      20.3     20.4 TRUE    holiday 
 4 2012-01-01 03:00:00  7256.      19.8     20.1 TRUE    holiday 
 5 2012-01-01 04:00:00  6793.      19.0     19.1 TRUE    holiday 
 6 2012-01-01 05:00:00  6636.      18.7     18.8 TRUE    holiday 
 7 2012-01-01 06:00:00  6548.      18.7     18.8 TRUE    holiday 
 8 2012-01-01 07:00:00  6865.      19.6     20.1 TRUE    holiday 
 9 2012-01-01 08:00:00  7300.      21.8     22.6 TRUE    holiday 
10 2012-01-01 09:00:00  8002.      24.6     25.2 TRUE    holiday 
# … with 26,294 more rows
p <- autoplot(vic_elec_hour, demand, color = "blue")

ggplotly(p, dynamicTicks = TRUE) |> rangeslider()
vic_elec_hour |> 
  ggplot(aes(x = mean_temp, y = demand, 
             color = holiday)) +
  geom_point(alpha = 0.5)

vic_elec_hour |> 
  ggplot(aes(x = mean_temp, y = demand, 
             color = day_type)) +
  geom_point(alpha = 0.5)

(train <- filter(vic_elec_hour, year(time) %in% 2012:2013))
# A tsibble: 17,544 x 6 [1h] <Australia/Melbourne>
   time                demand mean_temp max_temp holiday day_type
   <dttm>               <dbl>     <dbl>    <dbl> <lgl>   <chr>   
 1 2012-01-01 00:00:00  8646.      21.2     21.4 TRUE    holiday 
 2 2012-01-01 01:00:00  7927.      20.6     20.7 TRUE    holiday 
 3 2012-01-01 02:00:00  7902.      20.3     20.4 TRUE    holiday 
 4 2012-01-01 03:00:00  7256.      19.8     20.1 TRUE    holiday 
 5 2012-01-01 04:00:00  6793.      19.0     19.1 TRUE    holiday 
 6 2012-01-01 05:00:00  6636.      18.7     18.8 TRUE    holiday 
 7 2012-01-01 06:00:00  6548.      18.7     18.8 TRUE    holiday 
 8 2012-01-01 07:00:00  6865.      19.6     20.1 TRUE    holiday 
 9 2012-01-01 08:00:00  7300.      21.8     22.6 TRUE    holiday 
10 2012-01-01 09:00:00  8002.      24.6     25.2 TRUE    holiday 
# … with 17,534 more rows
# maximos:
# year = 4380
# week = 84
# day = 12

fit <- train |>  
  model(
    Fourier = TSLM(log(demand) ~ fourier(period = "year", K = 438) + 
                      fourier(period = "week", K = 42) + 
                      fourier(period = "day", K = 6)))
glance(fit)
# A tibble: 1 × 15
  .model  r_squa…¹ adj_r…²  sigma2 stati…³ p_value    df log_lik     AIC    AICc
  <chr>      <dbl>   <dbl>   <dbl>   <dbl>   <dbl> <int>   <dbl>   <dbl>   <dbl>
1 Fourier    0.876   0.869 0.00444    122.       0   961  23128. -94120. -94008.
# … with 5 more variables: BIC <dbl>, CV <dbl>, deviance <dbl>,
#   df.residual <int>, rank <int>, and abbreviated variable names ¹​r_squared,
#   ²​adj_r_squared, ³​statistic
fc <- forecast(fit, h = "1 year")
Warning: There was 1 warning in `mutate()`.
ℹ In argument: `Fourier = (function (object, ...) ...`.
Caused by warning:
! prediction from a rank-deficient fit may be misleading
fc
# A fable: 8,766 x 4 [1h] <Australia/Melbourne>
# Key:     .model [1]
   .model  time                           demand  .mean
   <chr>   <dttm>                         <dist>  <dbl>
 1 Fourier 2014-01-01 00:00:00   t(N(9, 0.0047))  8399.
 2 Fourier 2014-01-01 01:00:00   t(N(9, 0.0047))  8329.
 3 Fourier 2014-01-01 02:00:00   t(N(9, 0.0047))  8063.
 4 Fourier 2014-01-01 03:00:00 t(N(8.9, 0.0047))  7697.
 5 Fourier 2014-01-01 04:00:00 t(N(8.9, 0.0047))  7484.
 6 Fourier 2014-01-01 05:00:00   t(N(9, 0.0047))  7727.
 7 Fourier 2014-01-01 06:00:00 t(N(9.1, 0.0047))  8572.
 8 Fourier 2014-01-01 07:00:00 t(N(9.2, 0.0047))  9740.
 9 Fourier 2014-01-01 08:00:00 t(N(9.3, 0.0047)) 10537.
10 Fourier 2014-01-01 09:00:00 t(N(9.3, 0.0047)) 10569.
# … with 8,756 more rows
ggplot() +
  geom_line(vic_elec_hour, mapping = aes(x = time, y = demand)) + 
  geom_line(fc, mapping = aes(x = time, y = .mean), color = "blue", alpha = 0.8)